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INTRODUCTION 


Continuing  research  into  hybrid  stress  finite  element  formulations  has  been  directed  towards  minimizing  the 
computational  cost  associated  with  computing  element  stiffness  coefficients  in  order  to  make  hybrid  elements 
competitive  with  purely  displacement-based  formulations.  Much  emphasis  has  been  placed  on  reducing  the 
computational  cost  of  inverting  the  complementary  energy  or  flexibility  matrix  inherent  to  the  hybrid  stress 
method  and  has  led  to  insightful  yet  elaborate  procedures  and  approximations  directed  towards  this  aim^’^’^. 
Recently,  a  general  procedure  has  been  developed  by  the  author  for  simplifying  the  stiffness  definition  by 
making  full  use  of  the  freedom  in  selecting  and  manipulating  assumed  stress  fields  which  include  the  gener¬ 
ation  of  a  weighted  orthonormalized  stress  mode  basis  In  this  approach,  the  flexibility  matrix  is  formally 
eliminated  and  the  resulting  expression  for  the  element  stiffness  matrix  is  reduced  to  the  integration  of  a 
single  constituent  matrix  which  can  be  accomplished  in  closed-form.  Such  an  evaluation  results  in  algebraic 
equations  which  can  appeeur  cumbersome;  however,  these  statements  replace  the  formation  and  inversion  of 
the  flexibility  matrix  and,  more  importantly,  the  costly  numerical  quadrature  of  large-order  matrix  products. 
In  contrast,  displacement-based  continuum  elements  under  general  distortion  possess  an  integrand  involving 
rational  functions  due  to  the  presence  of  the  Jacobian  determinant  in  the  denominator  which  precludes  a 
simple  explicit  evaluation  of  stiffness  coefficients  in  algebraic  form.  Such  a  represention  of  basic  strain  energy 
characteristics  by  rational  functions  compounded  by  possible  nonlinear  variations  in  material  properties  make 
the  inherent  polynomial  approximation  of  Gaussian  quadrature  schemes  significant  in  regards  to  the  rate  of 
solution  convergence.  The  present  method  is  completely  generic  and  potentially  avoids  all  approximations 
made  to  the  formal  variational  definition  of  hybrid  element  stiffness  matrices.  The  developed  technique  is 
utilized  in  two  different  element  derivations.  First,  a  closed- form  set  of  expressions  are  developed  for  the 
stiffness  coefficients  in  the  Pian-Tong  hexahedral  element  ®  which  is  a  robust  8-node  solid  continuum  element 
and  presents  a  formidable  degree  of  difficulty  in  deriving  an  explicit  formulation.  Secondly,  the  method 
is  extended  to  derive  explicit  expressions  for  element  stiffness  matrices  incorporating  nonconst2int  material 
properties  for  nonlinear-elastic  problems  which,  within  the  framework  of  finite  element  solution  methods, 
usually  require  computationally  intensive,; iterative  solution  procedures.  As  an  initial  effort  towards  that 
aim,  an  explicit  formulation  is  developed  for  the  4-node  Pian-Sumihara  hybrid  quadrilateral  element®.  The 
resulting  closed-form  derivations  demonstrate  a  substantial  reduction  in  computational  cost  over  purely  nu¬ 
merical  treatments  in  generating  element  matrices. 

VARIATIONAL  BASIS  OF  THE  HYBRID  MODEL 

The  form  of  the  Hellinger-Reissner  energy  functional  utilized  in  References  [5]  and  [6]  is  given  by 

[(-l/2)o‘^S<r  +  <r'^{Luq)  -  (r)'^ u\]dv  (1) 

where  <r  is  the  assumed  stress  field,  S  is  the  material  compliance  matrix,  and  ua  are  the  assumed  compat¬ 
ible  and  incompatible  displacement  fields  and  L  is  the  differential  operator  relating  strains  to  displacements. 
The  assumed  stresses  may  be  represented  by 

a  =  P0  (2) 

where  P  is  a  matrix  of  polynomial  terms  and  /3  is  a  vector  of  undetermined  expansion  coefficients.  The 
displacement  field  is  assumed  over  the  element  domain  as 

u  =  -huA  =  Nq-hMA  (3) 

where  N  and  M  are  compatible  and  incompatible  displacement  shape  functions,  respectively,  q  are  nodal 
displacements,  and  A  are  Lagrange  multipliers  which  variationally  enforce  the  field  equilibrium  conditions. 
These  variational  constraints  are  applied  a  priori  to  the  assumed  stress  modes  which  condense  the  influence 
of  the  incompatible  modes  into  the  element  formulation.  The  resulting  expression  for  the  element  stiffness 
matrix  is  given  by 

K  =  G^H-^G  (4) 

where 

H  =  /  P'^SPdv  (5) 

Jv 


1 


-L 


(6) 


in  which  B  is  the  strain-displacement  matrix.  In  the  above  formulation,  the  a  priori  condensation  of  the 
incompatible  displacement  modes  results  in  an  H  matrix  which  is  fully  populated  and  has  hitherto  required 
the  inversion  of  a  full  matrix  of  order  dim{l3). 

COMPUTATION  OF  EXPLICIT  STIFFNESS  MATRICES 

The  procedure  developed  in  Reference  [4]  for  simplifying  the  expressions  involved  in  (4)  is  outlined  below. 
The  method  utilizes  a  sequence  of  permissible  transformations  of  assumed  stress  fields  to  simplify  the  stiffness 
definition.  The  assumed  stresses  are  first  transformed  through  the  introduction  of  a  symmetric  ‘distributing’ 
matrix  as 

P=IP  =  (DD-^)P  =  D(D-‘P)  =  DP  (7) 

where  the  distributing  matrix  is  defined  as 

D  =  S-i/2  (8) 

The  inverse  square  root  of  the  compliance  matrix  is  obtained  through  a  standard  spectral  decomposition^.  For 
illustration,  the  decomposition  of  2-D  and  3--D  orthotropic  compliance  matrices  are  presented  in  Appendix 
I.  Substitution  of  (7)  into  (5)  yields 


H  =  y  j  j  [|J|P^D^SDP](i|d>7dC 


(9) 


where,  from  the  definition  of  D  and  the  symmetry  of  both  S  and  D,  we  obtain 

D^SD  =  =  SS”^  =  I 

and  the  flexibility  matrix  reduces  to 


K  =  J'j'j'[\3\P'^P]d^dridC 


(10) 


A  second  field  transformation  uses  equation  (10)  to  define  a  weighted  inner  product  for  use  in  a  Gram- 
Schmidt  procedure  to  generate  an  orthonormal  spanning  set  of  stress  modes,  P'*',  which  are  a  special  linear 
combination  of  the  modes  present  in  P.  The  weighted  inner  product  is  therefore  defined  as 


P.-.P;  >=  fj^y3\PlPiHdr,dC  =  Sij 


(11) 


where  is  the  Kronecker  delta  function.  The  linear  combination  yielding  a  sequence  of  orthogonal  stress 
modes  is  given  by 

^  P,-  i  =  1 

'  p<  -  <  Pj*. p.-  >  p)  » >  1 

.  i=i 

which  are  normalized  to  form  basis  vectors,  P*,  as 


pr  =<  Vi,  Vi  y. 

Substitution  of  P’  into  equation  (10)  yields  by  definition 

H  =  y' [|J|P*^P*]ded;?dC  =  I 

The  new  basis  for  the  element  stress  field  is  now  given  by 

P  =  DP* 


(13) 


(14) 


(15) 
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and  the  expression  for  the  element  stiffness  matrix  reduces  to 


K  =  G^G 


Separating  out  the  Jacobian  determinant  from  the  isoparametric  strains  as 

B= 

and  substituting  (15)  and  (17)  into  (6)  yields  the  G  matrix  definition  as 


G  =  y  jP*^DB*](iedr/dC 


(16) 


(17) 


(18) 


The  absence  of  the  Jcicobian  determinant  in  the  denominator  permits  a  direct  derivation  of  algebraic  expres¬ 
sions  for  the  G  matrix  coefficients.  The  explicit  form  of  the  element  stiffness  matrix  is  then  obtained  from 
equation  (16). 


EXPLICIT  PIAN-TONG  HEXAHEDRAL  ELEMENT 


The  configuration  of  the  8-node  Pian-Tong  hex2Lhedral  element  is  depicted  in  Figure  1.  The  compatible 
displacement  functions  Uq  are  given  by 


(19) 


As  detailed  in  Reference  [5],  incompatible  displacement  modes  are  introduced  to  complete  the  cubic  order 
of  the  assumed  isoparametric  displacement  field.  These  modes  are  condensed  a  priori  into  the  element  for¬ 
mulation  through  constr Clint  conditions  on  the  assumed  stress  modes. 


The  stress  field  utilized  in  the  Pian-Tong  element  is  assumed  in  natural  coordinates  cind  transformed  to 
physical  coordinates  using  Jacobians  computed  at  the  element  centroid  as 

=  (20) 

Performing  the  initial  transformation  of  stresses  given  by  (7)  results  in  the  physical  or  Cartesian  stress  field 
given  by 

(T  =  [cTxiO’yyO’z^Tyz^TzxiTxy]  (2^) 
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where 


and 


Pk  = 
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(22) 


(23) 


The  coefficients  are  obtained  from  the  product  of  the  inverse  distributing  matrix  and  a  transformation 
matrix  relating  tensorial  stresses  to  physical  components.  This  relationship  is  given  by 


E  =  D^^T 


or 


■  ^11 

di2 

<^13 

a? 

al 

«3 

2aia2 

20203 

20103 

^21 

^22 

^23 

6? 

6| 

ii 

26162 

26263 

26163 

[e.;]  = 

dai 

^32 

^33  ^ 

c? 
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C§ 
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aiCi 
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The  isoparametric  mapping  between  local 

physical  and  natural 

coordinates  is 

given  by 

X  —  +  a2f}  + 

y  =  +  ^>21  +  H"  (24) 

+  C2V  +  caC  +  C4^V  +  cs^C  +  ceVC  +  <^7^vC 

where 
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Orthonormalizing  the  stress  modes  according  to  (12)  yields  the  final  transformed  stress  field  composed  of 
constant  modes  given  by 


K  =  niPc  (25) 

and  higher-order  modes  given  by 

“  [P7jP8»P9)Pl0>Pll)Pl2Jpia>Pl4jPl5»  PiSjPiTjPis]  (26) 
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The  columns  of  the  higher^rder  stress  modes  with  i  =  1, 2, 3, 6  and  j  =  1,2,3  are  given  by 
{Pij+e}  =  «i+6{r<,2i-i^  +  n',2i} 

{Pr,;+9}  =  ”i+9{'*«,3i+4»7  +  »’i,3;+5^  +  n.Sj+e} 

{Pi,i+ 12}  =  "j+12{»*«,4j  +  12C  +  »’»,4i+13»?  +  '’»,4i  +  144  +  ’’«,4i  +  15} 

{p.*!i6}  =  ni6{ri,28^»7  +  r<,29C  +  n,30»7  +  +  r,-,32} 

{Pi,!?}  nni^i, 33T]C  +  ''•,34^»7  +  Ti.asC  +  **»,36»7  +  ^i,37^  +  »*,-,38} 

{Pi.is}  =  ni8{»*«,39€C  +  ^i,40VC  +  '■i,41^»?  +  »*«,42C  +  ri^^T]  +  ^-,44^  +  r,-,45} 

where  the  expressions  for  n,*  and  r,  j  are  presented  in  Appendix  II. 

The  integration  of  (18)  yields  the  following  expressions  for  elements  of  the  G  matrix  in  which  j  =  1, 2, 3, 8. 
Specific  components  utilize  an  additional  index,  i,  which  assumes  the  values  i  =  1,2, 3.  The  components  are 
given  by 

9i,3j-2  =  di,Ti(l,j)  gi,3j-i  =  d2,Ti(2,j)  gi^3j  —  d3iri(3,j) 

^4,3i-2  =  0  94,3j--l  =  d44ri(3,j)  94, 3j  =  d44ri(2,Ji‘) 

95,3j^2  =  dss^i{ZJ)  =  0  95, 3j  =  C^55ri(l,;) 

^6,3i-2  =  dee^lC^J)  96,3j-l  =  <^66ri(l,i)  96, 3j  =  0 

9i4-6,3j--2  =  <^iir'2(li4,  i,  j)  +  di2r2(2,4,  i,i)  4-  di3r2(3,4,  i,  j)  +  d55r2(5, 6,  i,  j)  +  d66r2(6,5,  i^j) 

9i+6,3j^l  =  <^2ir2(l,5,  i,i)  +  d22r2(2,5,  i,j)  +  d23r2(3,5,  l,i)  +  d44r2(4,  6,i,  j)  +  d66r2(6,4,  1,  j) 

flf»+6,3i  =  d3ir2(l,6,i,  j)  4- d32r2(2,6,  i,  j)  4- d33r2(3,6,i,  j)  +  d44r2(4, 5,  i,  j)  4- d55r2(5,4, 2,i) 

9i-{~9,3j-2  —  diir3(l,7,  z,  j)  4-  di2r3(2, 7,  z,  j)  4-  di3r3(3,7,  z,  j)  4-  (issTaCS,  9,  z,  j)  4-  d66r3(6,8,  z,  j)  (27) 

5^1+9, 3i-l  =  <^2ir3(l,8,  Z,  j)  4- d22r3(2,8,  Z,  j)  4- d23r3(3,8,  Z,  j)  4‘d44r3(4,9,Z,  j)  4- d66r3(6,7,  z,  j) 

^•+9,3i  =  d3ir3(l,9,  z,  j)  4-  d32r3(2, 9,  z,  j)  4-  d33r3(3,9,  z,  j)  4“  d44r3(4, 8,  z,  j)  4-  d55r3(5,7,  z,  j) 

9i4-l2,3j-2  —  diir4(l,  10,  Z,  j)  4-  di2r4(2,  10,  i,j)  4“  di3r4(3,  10,  Z,i)  +  d55r4(5,  12,  Z,^')  4“  ^661^4(6,  11,  i^j) 

^t+i2,3i~i  =  d2ir4(l,  11,Z,  j)  4-  d22r4(2,  ll,z,  j)  4-  d23r4(3, 11,  ij)  4-  d44r4(4, 12,  ij)  4-  10,z,  j) 

9i+12,3j  =  d3ir4(l,  12,  Z,  j)  4-  d32r4(2,  12, 1,^)  4-  d33r4(3,  12,  iJ)  4-  d44r4(4,  ll,z,  j)  +  d55r4(5,  10,.z,  j) 

flfi6,3j-2  =  ^^iir5(l,  13, i)  4-  di2r5(2, 13,  j)  4“  di3r5(3, 13,  j)  4-  d55r5(5, 15,  j)  4-  deeVsi^,  14,  j) 

^16,3;-l  =  <^2ir5(l,  14, i)  4-  d22r5(2,  14,  j)  4"  d23r5(3,  14,  j)  4"  d44r5(4,  15,  j)  4-  d66r5(6,  13, i) 

</l6,3i  =  ^^3ir5(l,  15, i)  4-  d32r5(2,  15,  j)  4-  d33r5(3,  15,  j)  4-  d44r5(4, 14,  j)  4*  d55r5(5,  13,  j) 

^i7,3j-2  =  diir6(l,  16,  j)  +  di2r6(2, 16,  j)  4-  di3r6(3, 16,^')  +  d55r6(5, 18,  j)  4-  d66r6(6, 17,  j) 

^17,3;-!  =  d2ir6(l,  17,  j)  4*  d22r6(2,  17,  j)  4“  d23r6(3,  17,  j)  4“  d44r6(4,  18, i)  4-  d66r6(6,  16,  j) 

^17,3;  =  <^3ir6(l,  18, i)  4*  d32r6(2,  18,  j)  +  d33r6(3,  18,  j)  4"  d44r6(4,  17,  j)  4*  d55r6(5,  16,  j) 

5^18, 3i-2  =  <iiir7(l,  19, i)  4-  di2r7(2, 19,  j)  4-  di3r7(3, 19,  j)  +  d55r7(5, 21,  j)  4-  d66r7(6,20,  j) 

9is,3j-i  =  d2ir7(l,20,  j)  4*  d22r7(2,20,  j)  4-  d23r7(3,20,  j)  4-  d44r7(4, 21,  j)  4-  dee^ii^^  19,  j) 

9is,3j  =  d3ir7(l,21,  j)  4*  d32r7(2,21,  j)  4-  d33r7(3, 21,  j)  +  d44r7(4, 20,  j)  4-  d55r7(5, 19,  j) 

where  dij  are  elements  of  the  distributing  matrix  and  the  functions  F,-  are  given  by 

ri(m,i)  =  ni<^j,m/3 

r2(fc,  m,  Z,  j)  =  ni^e{n,2i-l<l>j,m  +  rjb.2i<^i,m-3)/9 

r3(A:,m,Z,  j)  =  Tli^^{rk^3i4.4<l>j,m  4- rjt,3j4.5<^j,m-3  4- rfc^3,*4.6<^j,m-6)/9 

T4(k,  7n,iyj)  =  ni^l2{^k,4i+l2<l>j,m  4-  rjb,4,-+i3<^j,ni-3  4*  rjb,4t+14<^i,m-6  4-  rk,4i4'lB<i>j,m^9)/^ 

Ts{ky  m,  j)  =  nisirk,2S<l>j,m  4-  rk,29(l>j,m-3  +  rk,3Q<l>j,m^6  4-  rjfc,31<^i,m~9  4-  rjb,32<^j,m-.12)/27 

r6(fc,  rriyj)  =  ni7{rk,33<l>j,m  4-  ,34<pj ,m--3  4*  rk,35<l>j,m^6  4-  rk,3e</>j,m-9  4-  rk,37<l>j,m-12  4-  ri,38<^j,m-15)/27 

r7(fc,  rriyj)  =  nis(rk,39<t>j,m  4-  rk,40<l>j,m^3  4-  rk,4l<f>j,m^6  4*  rfc,42<^j,m-9  4-  Z*fc,43<^i,m-124- 

^fc,44<^i,m-15  +  ^jfc,45<^i,m-18)/27 
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in  which  geometric  constants 


1  r»  1 


6lj  =  biCj’-bjCi  6?-  =  CiOj  —  CjGi 

and  values  for  with  n  =  1,2,3,  ..,,8  are  given  by 


CL^bj  CL 2  b^ 


n 

^3 

Z4 

-n 

Z? 

1 

-1 

1 

-1 

1 

-1 

-1 

1 

2 

1 

-1 

1 

-1 

-1 

-1 

1 

3 

1 

-1 

-1 

1 

-1 

1 

-1 

4 

-1 

1 

1 

-1 

-1 

1 

-1 

5 

-1 

-1 

-1 

-1 

1 

1 

1 

6 

1 

1 

1 

1 

1 

1 

1 

7 

1 

1 

-1 

-1 

1 

-1 

-1 

8 

-1 

-1 

1 

1 

1 

-1 

-1 

(28) 


NUMERICAL  STUDIES  OF  THE  EXPLICIT  PIAN-TONG  ELEMENT 

The  following  tables  present  computer  run-times  comparing  the  explicit  and  numerical  generation  of  stiffness 
matrices  in  the  Pian-Tong  hexahedral  element.  In  the  numerical  evaluation,  under  general  element  distor¬ 
tion,  quartic  terms  involving  the  natural  coordinates  appear  requiring  a  3^‘^-order  Gaussian  quadrature  rule 
for  exact  evaluation  as  is  obtained  in  the  explicit  formulation.  However,  the  under  integration  utilizing  a  2"^- 
order  rule  has  been  found  adequate  to  compute  element  stiffness  coefficients  without  introducing  spurious 
kinematic  modes.  The  influence  of  the  quartic  terms  have  a  negligible  influence  and  a  2”‘^-order  rule  is  thus 
considered  standard  for  the  8-node  solid  hybrid  element  and  is  used  herein  for  all  numerical  integrations. 
No  optimization  was  attempted  in  performing  the  computer  implementations  of  the  elements  in  terms  of 
code  preparation  or  CPU  processing  options  such  as  vectorization  or  concurrency.  The  codes  were  run  on  a 
Hewlett  Packard  Apollo  400  series  workstation  in  a  Unix  environment.  The  standard  UNIX  profiler  Gprof 
was  used  to  characterize  the  time  spent  in  performing  various  operations.  In  generating  the  computational 
profiles,  1000  element  stiffness  matrices  were  processed.  A  description  of  the  procedures  used  in  the  codes 
which  account  for  at  least  97%  of  the  computational  cost  are  presented  in  Table  1.  The  designation  ^main’ 
combines  the  operations  within  the  main  program  together  with  various  subroutines  which  contribute  in¬ 
significant  computational  cost.  Table  2  details  the  computational  profile  for  the  Pian-Tong  element  evaluated 
explicitly  and  numerically. 


Table  1:  Subroutine  procedures 


Name 

Description 

Name 

Description 

mxmul 

mxadd 

invers 

matrix  multiplication 
matrix  addition 
matrix  inversion 

gmatrx 

orthop 

main 

explicit  computation  of  G  matrix 
computation  of  orthonormal  stress  modes 
main  program  -h  minor  subroutines 
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Table  2:  Computational  profiles  of  the  Pian-Tong  (PT)  element 


PT- Explicit 

PT-Numerical 

% 

self 

cumulative 

procedure 

% 

self 

cumulative 

procedure 

time 

seconds 

seconds 

name 

time 

seconds 

seconds 

name 

70.5 

39.99 

39.99 

mxmul 

83.5 

341.62 

341.62 

mxmul 

18.6 

10.56 

50.55 

gmatrx 

7.3 

29.90 

371.52 

mxadd 

8.9 

5.03 

55.58 

orthop 

6.9 

28.43 

399.95 

invers 

2.0 

1.12 

56.70 

main 

2.3 

9.42 

409.37 

main 

The  computational  profiles  quantify  the  different  characteristics  of  the  explicit  and  numerical  versions  of 
the  Pian-Tong  element.  As  shown,  the  computational  cost  of  both  the  explicit  and  numerical  versions  is 
mostly  spent  in  matrix  multiplication  and  addition  consuming  70.5  %  and  90.8  %  of  the  processing  time, 
respectively.  Comparing  total  computer  run-times,  the  explicit  version  requires  only  13.9  %  of  the  compu¬ 
tational  cost  as  that  consumed  in  the  underintegrated  numerical  evaluation.  Not  shown  in  the  tables,  the 
exact  integration  of  the  quartic  terms  in  the  numerical  version  using  3^^-order  Gaussian  quadrature  results 
in  a  run-time  of  1167.4  seconds.  Compared  to  this  the  explicit  formulation  requires  only  4.9  %  of  the  com¬ 
putational  expense. 


NONLINEAR  ELEMENT  STIFFNESS  DEFINITION 


The  extension  to  nonconstant  material  properties  requires  that  the  compliance  matrix  be  interpolated  over 
the  element  domain.  For  illustration,  the  Pian-Sumihara  quadrilateral  element  is  selected^.  For  the  plane  ele¬ 
ment  under  study  a  bilinear  isoparametric  field  is  sufficient  for  interpolating  compliance  properties;  however, 
the  derived  D  matrix  is  a  nonlinear  function  of  S  thus  requiring  a  higher  order  interpolation.  A  quadratic 
isoparametric  field  is  adopted  in  which  the  material  matrices  are  computed  at  the  stress  recovery  or  material 
evaluation  points  depicted  in  Figure  2.  Thus,  at  arbitrary  points  over  the  element  domain  D  =  and 

the  interpolation  results  in  the  relationship  S  ==  S  being  strictly  valid  at  the  recovery  points  and  approximate 
elsewhere.  The  accuracy  of  this  approximation  is  discussed  later  in  Rem2ULk  I. 


Figure  2.  Evaluation  points  used  for  material  interpolation. 


The  interpolation  is  given  by 


8 

i=l 


(29) 


where  /?,-  represents  the  material  D  and  matrices  computed  at  the  stress  recovery  point  and  Ni  are 
the  associated  quadratic-order  isoparametric  shape  functions.  Interpolation  of  the  material  matrices  over 
the  element  domain  may  be  expressed  as 


v)  —  +  A2^  +  Aar;  4-  A4^^  4-  Ast;^  4-  Ae^r;  4-  Ay^^r;  4-  Ag^r;^  (30) 
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where  the  constant  matrices,  A,*,  are  given  by 


f  1 

■  _1  _1  _1  _1  2  2  2  2  ■ 

'  /?!  " 

A2 

2  -2 

(I2 

M 

-2  2 

1?3 

A4 

_  1 

1-1  1-1 

Ha 

As 

^  ~  4 

1111-2  -2 

fl. 

Ae 

-1-111  2  -2 

At 

-1-1112  -2 

flj 

.  As  , 

-111-1  -2  2 

.  ■^8  , 

(31) 


Following  the  procedure  outlined  above  in  equations  (7)  through  (18),  substitution  of  (7)  into  (5)  based  on 
the  distribution  represented  by  equation  (30)  yields  the  flexibility  matrix  as 


H  =  y  j  [|J|P^D'^SDP]d^d7; 


where,  from  the  definition  of  D  and  the  symmetry  of  both  S  and  D  we  obtain 

D^SD  =  S-l/2gg-l/2  _  gg-l  ^  j 


(32) 


(33) 


The  approximation  rapidly  approeiches  an  identity  as  the  interpolation  order  of  D  increases  and  as  the  varia^ 
tion  of  material  properties  over  the  element  diminishes.  The  flexibility  matrix  is  thereby  closely  approximated 

by  .  ^  ^ 

/  /  [IJIP^Pjd^dr,  (34) 


A  second  field  transformation  uses  equation  (34)  to  define  a  weighted  inner  product  for  use  in  a  Gram- 
Schmidt  procedure  to  generate  an  orthonormal  spanning  set  of  stress  modes,  P*.  This  modal  set  is  formed 
as  a  special  linear  combination  of  the  modes  present  in  P  and  constitutes  an  alternative  basis  set  which 
spans  the  original  stress  space.  The  weighted  inner  product  is  therefore  defined  as 


<  Pi,Pj  >=  j'j\3\PjPj]d^dr,  =  6, 


(35) 


where  Sij  is  the  Kronecker  delta  function.  The  linear  combination  yielding  a  sequence  of  orthogonal  stress 
modes  is  given  by 

P^  i  =  1 

i  —  1 

Pi  -  <  p;, Pi  >  p;  i  >  1 

i=i 

which  are  normalized  to  form  basis  vectors,  P*,  as 


Vi  = 


p;=<  Vi,Vi  > 


-1/2 


Vi 


(37) 


In  the  definition  of  P,  the  D“^  matrix  is  interpolated  according  to  (30)  and  carried  into  the  subsequent 
orthonormalization  of  assumed  stress,  modes.  For  the  present  study,  however,  a  simplification  is  adopted 
wherein  the  D  matrix  is  interpolated  according  to  equation  (30)  while  its  inverse  is  obtained  from  an  area 
average  obtained  by  integrating  over  the  element  domain.  This  operation  is  not  required  but  is  utilized  in 
the  present  study  to  simplify  the  expressions  for  the  orthonormalized  stress  modes,  P*,  while  maintaining 
an  adequate  approximation  for  the  material  properties.  The  inverse  is  thus  defined  as 


and  the  relationship 


D(^,r?)D-'«I 


-.-1 


(38) 

(39) 
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is  satisfied  in  an  integral  sense  for  nonconstant  material  properties  and  pointwise  for  constant  properties 
over  the  element  domain.  Substitution  of  P*  into  equation  (34)  yields  by  definition 

H  =  [|J|P*^P*]cie</77  =  I  (40) 

Repeating  the  above  relations,  the  new  basis  for  the  element  stress  field  is  given  by 

P  =  DP"  (41) 

and  the  expression  for  the  element  stiffness  matrix  reduces  to 

K  =  G^G  (42) 

Separating  out  the  Jacobian  determinant  from  the  isoparametric  strains  as 

B  =  (43) 

and  substituting  (41)  and  (43)  into  (6)  the  G  matrix  definition  becomes 

G  =  J  J  [P*’’DB*K(f7;  (44) 

The  absence  of  the  Jax^obian  determinant  in  the  denominator  permits  a  direct  derivation  of  algebraic  ex¬ 
pressions  for  the  G  matrix  coefficients  which  incorporate  a  nonconstant  field  of  material  properties  over  the 
element  domain.  The  explicit  form  of  the  element  stiffness  matrix  is  then  obtained  from  equation  (42). 


NONLINEAR  PIAN-SUMIHARA  QUADRILATERAL  ELEMENT 


Figure  3.  Quadrilateral  element  configuration. 


The  configuration  of  the  Pian-Sumihara  element  eind  node  numbering  are  depicted  in  Figure  3.  The  dis¬ 
placement  functions  are  given  by 

As  detailed  in  Reference  [6],  stresses  are  defined  in  natural  or  tensorial  coordinates  and  incompatible  dis¬ 
placement  modes  are  introduced  to  complete  the  quadratic  order  of  the  assumed  isoparametric  displacement 
field.  These  modes  are  condensed  a  priori  into  the  element  formulation  through  constraint  conditions  on  the 
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assumed  stress  modes. 


The  tensorial  stress  field  is  transformed  to  physical  or  Cartesian  coordinates  using  Jacobians  computed 
at  the  element  centroid  as 

(46) 

Performing  the  initial  transformation  of  stresses  given  by  (7)  results  in  the  Cartesian  stress  field  given  by 


where 


<T  ~  \<Tx ,  (Ty ,  Txy^  — 

1  ci2l 

P  =  I  1  C2l^  C22T} 

1  0327/ 


(47) 


(48) 


and  the  coefficients  are  obtained  from  the  product  of  D“^  and  a  transformation  matrix  relating  tensorial 
stresses  to  Cartesian  components.  This  relationship  is  given  by 


[c,,]  =  D”^r  = 


dll  di2 
^21  ^22 


^33 


4 

«i 

H 

^2^2 

aibi 

The  geometric  parameters  a,-  and  are  obtained  from  the  mapping  between  physical  and  natural  coordinates 
given  by 

X  =  ao  +  ai<f  H- 027/ +  03^77 

y  =  60  "h  bi^  4-  627/  -h  bs^T] 

where 


(49) 


’  ao 

bo  ' 

1 

1 

1 

1 

di 

bi 

_  1 

-1 

1 

1 

-1 

d2 

b2 

"  4 

-1  - 

1 

1 

1 

. 

63 

1  - 

1 

1 

-1 

yi  * 

^2 

y2 

^3 

ys 

X4 

P4  . 

The  weighted  orthonormalized  stress  modes  are  obtained  as 


P"  = 


Pl3  0  0  P42^  +  P43  PsI^+pSs^ +P53 

0  P26  0  Phi  +  Pie  P54^-bP55e  +  P56 

0  0  P39  P48^  +  P49  Plr^  +  Pss^  +  ph 


(50) 


where  the  stress  mode  coefficients,  are  presented  in  Appendix  II. 

Given  the  following  constants  arising  from  the  regular  structure  of  the  strain  modes 


Cij  —  ^2^  ””  b\Z2  ^4j 

62;  =  63^-614  ^5; 

63;  =  622:3  -  63^  66; 


aiz{^a2z{ 

aiz^  —  a3Z^  where 

(I3Z2  ~  ^2^ 


and  the  general  form  for  each  stress  mode  given  by 


j 

4 

~T 

4 

1 

-1 

-1 

1 

2 

1 

-1 

-1 

3 

1 

1 

1 

4 

-1 

1 

-1 

{Pn77  +  p,-2^-fpi3  ] 
Pt477  +  P,*5^  H-Pi6  > 
Pt77/+Pi8^ +  Pi9  J 


(51) 


(52) 


the  integration  of  (44)  yields  explicit  expressions  for  the  components  of  the  G  matrix.  The  expressions 
are  based  on  the  most  general  form  of  orthonormalized  stress  modes  in  which  most  terms  are  zero  for 
the  simpler  modes  and  the  constants  {dij)k  are  elements  of  the  D  matrix  computed  at  the  evaluation 
point.  In  addition,  for  constant  material  properties,  if  all  (d,;)jt  terms  are  discarded  for  ib  >  1  the  resulting 
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expressions  reduce  to  an  explicit  form  of  the  linear-elastic  Pian-Sumihara  element.  The  components  gij  with 
z  =  1, 2, 3, 5  and  j  =  1, 2, 3, 4  are  given  by 

=  (Pa{(^ij[(^n)7  +  3(dii)3]  +  e3j[(dii)4  -f  3(dii)i  +  9(dii)5/5]  +  (^11)6^2;}  + 

Pi2{(^ij[(^ii)8  +  3(dn)2]  +  e2j[(dii)5  +  3(dii)i  +  9(dii)4/5]  +  (^11)663^}  + 

Pi3{(^2j[(^ii)8  +  3((iii)2]  +  e3j[(dii)7  +  3(dii)3]  +  3eij[(dii)5  +  (dn)4  +  3(dii)i]}  + 
Pi4{(^lj[(^12)7  +  3(di2)3]  +  €3j[(dl2)4  +  3(di2)l  +  9(di2)5/5]  +  e2j(dl2)6}  + 

P?5{(^li[(^12)8  +  3(^12)21  +  e2i[(dl2)5  +  3(dl2)l  +  9(di2)4/5]  -f  e3j(di2)6}  + 

Pi6{(^3i[(^^12)7  +  3(di2)3]  +  ^2j[{dl2)s  +  3(dl2)2]  +  3ei;[(di2)5  +  (^12)4  +  3((ii2)l]}  + 
Pi7{(^4j[(<^33)7  +  3(6/33)3]  +  e6j[(d33)4  +  3(^33):  +  9(d33)5/5]  +  65^(6/33)6}  + 

Pi8{(^4j[(<^33)8  +  3(6/33)2]  +  e5j[(d33)5  +  3(6/33)1  +  9(d33)4/5]  +  e6y(6/33)6}  + 

Pi9{(^5j  [(6/33)8  +  3(6/33)2]  H-  eej [(6/33)7  +  3(6/33)3]  +  3e4;  [(6/33)5  4-  (6/33)4  +  3(6/33)i]})/9 

(53) 

9i,2j  =  (Pil{(^4j  [(6/12)7  +  3(6/12)3]  +  66i  [(6/12)4  +  3(6/12)1  +  9(6/12)5/5]  +  (6/12)665;  }  + 

Pi2  {(^4;  [(6/12)8  +  3(6/12)2]  +  65j  [(6/12)5  +  3(6/i2)l  +  9(6/12)4/5]  +  (6/12)666;}  + 

P<3  {(^5;  [(6/12)8  +  3(6/12)2]  +  66;[(6/i2)7  +  3(6/12)3]  +  3e4j[(6/i2)5  +  (6/12)4  +  3(6ii2)l]}  + 

P«*4  {(^4;  [(6/22)7  +  3(6/22)3]  H-  66;[(6/22)4  4*  3(6/22)1  +  9(6/22)5/5]  +  65^(6/22)6}  + 

Pi5  {(^4;  [(6/22)8  +  3(6/22)2]  +  65;  [(6/22)5  +  3(6/22)1  4  9(6/22)4/5]  +  66;  (6/22)$}  + 

Pi*6 {(^6;  [(6/22)7  +  3(6/22)3]  4-  e5j[(6/22)8  +  3(6/22)2]  4*  3e4j[( 6/22)5  +  (^22)4  +  3(6/22)1]}  + 
Pi7{(^lj[(^33)7  4*  3(6/33)3]  4  63j  [(6/33)4  +  3(6/33)1  4  9(6/33)5/6]  +  62;  (6/33)5}  4 

Pr8{(^lj[(6/33)8  4*  3(6/33)2]  +  e2;[(6/33)5  -f  3(6/33)1  -f  9(6/33)4/6]  +  63^(6/33)5}  4 

Pi's  {(^2;  [(6/33)8  4“  3(6/33)2]  +  e3j[(6/33)7  -f*  3(6/33)3]  +  36ij[(6/33)5  (6/33)4  +  H^33)l]})/^ 

NUMERICAL  STUDIES  OF  THE  EXPLICIT  PIAN-SUMIHARA  ELEMENT 

Computer  codes  were  generated  to  assess  element  computational  characteristics.  A  description  of  the  proce¬ 
dures  used  in  the  codes  are  presented  in  Table  3.  The  designation  ‘main’  combines  the  operations  within  the 
main  program  together  with  various  minor  subroutines  which  contribute  insignificant  computational  cost.  No 
optimization  was  attempted  in  terms  of  code  preparation  or  CPU  processing  options  such  as  vectorization  or 
concurrency.  The  codes  were  run  on  a  Hewlett  Packsird  Apollo  400  series  workstation  in  a  Unix  environment. 
The  standard  Unix  profiler  Gprof  was  used  to  characterize  the  time  spent  in  performing  various  operations. 
Table  4  presents  computational  profiles  and  computer  run-times  comparing  the  explicit  and  numerical  gener¬ 
ation  of  stiffness  matrices  in  the  nonlinear  Pian-Sumihara  quadrilaterzJ  element.  An  additional  comparison  is 
made  to  a  4-node  displacement-based  element  incorporating  incompatible  displacement  modes  which  is  pre¬ 
sented  in  Table  5.  The  incompatible  modes  are  based  on  quadratic  functions  and  modified  using  a  technique 
presented  in  Reference  [8]  to  identically  satisfy  the  strong  form  of  the  patch  test  for  incompatible  elements. 
A  2”‘^-order  Gaussian  quadrature  rule  was  used  for  all  numerical  evaluations  in  generating  computational 
profiles.  In  generating  the  computational  profiles,  10, 000  element  stiffness  matrices  were  processed. 


Table  3:  Subroutine  procedures 


Name 

Description 

Name 

Description 

mxmul 

mxadd 

invers 

state 

matrix  multiplication 
matrix  addition 
matrix  inversion 
static  condensation 

gmatrx 

orthop 

spectrl 

main 

explicit  computation  of  G  matrix 
computation  of  orthonormal  stress  modes 
spectra]  decomposition  of  S  matrix 
matrix  main  program  +  minor  subroutines 
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Table  4:  Computational  profiles  of  the  nonlinear  Pian-Sumihara  (PS)  element 


PS-Explicit 

PS-NumericaJ 

% 

self 

cumulative 

procedure 

% 

self 

cumulative 

procedure 

time 

seconds 

seconds 

name 

time 

seconds 

seconds 

name 

32.76 

16.55 

16.55 

spectrl 

75.46 

109.09 

109.09 

mxmul 

30.09 

15.20 

31.75 

mxmul 

9.26 

13.38 

122.47 

mxadd 

21.71 

10.97 

42.72 

gmatrx 

6.50 

9.39 

131.86 

in  vers 

3.17 

1.60 

44.32 

orthop 

8.78 

12.71 

144.57 

main 

12.27 

6.20 

50.52 

main 

Table  5:  Computational  profile  of  incompatible  displacement-based  element 


1  D-Based 

% 

self 

cumulative 

procedure 

time 

seconds 

seconds 

name 

71.6 

160.82 

160.82 

mxmul 

13.0 

29.18 

190.00 

mxadd 

11.0 

24.63 

214.63 

state 

4.4 

10.12 

224.75 

main 

The  computational  profiles  qucintify  the  different  characteristics  of  the  explicit  and  numerical  versions  of 
the  nonlinear  Pian-Sumihara  element.  In  the  explicit  version,  the  eight  spectral  decompositions  of  the  com¬ 
pliance  matrix  consume  the  greatest  amount  of  computational  cost  (32.76%)  while  matrix  multiplications 
constitute  most  of  the  computations  in  the  numerical  version  (75.46%).  The  computational  profile  of  the 
explicit  version  shows  that  the  operations  involved  in  forming  the  orthonormal  stress  modes  is  insignificant 
(3.17%)  while  the  formation  of  the  G  matrix  consumes  21.71%  of  the  cost.  The  final  evaluation  of  equation 
(24),  which  is  the  only  matrix  operation  performed  in  the  explicit  version,  constitutes  fully  30%  of  the  cost 
in  forming  the  element  stiffness  matrix.  Comparing  total  cost,  the  explicit  version  requires  only  34.95%  of 
the  processing  time  as  the  numerical  evaluation  and  only  22.48%  of  the  cost  required  in  the  incompatible 
displacement-based  formulation.  While  this  represents  a  significant  reduction  in  processing  cost,  application 
of  the  developed  methodology  may  be  expected  to  show  greater  reductions  in  3-D  and  higher-order  hybrid 
element  formulations.  This  issue  is  discussed  in  Remark  II. 

A  second  demonstration  is  made  to  show  the  accuracy  of  the  nonlinear  material  representation  in  the  ex¬ 
plicit  formulation.  Because  the  basic  computational  characteristics  have  been  shown  above,  a  cantilevered 
beam  under  plane  stress  is  solved  where  the  material  properties  are  not  a  function  of  the  stress /strain  state 
but  instead  vary  along  the  length  of  the  beam.  Such  a  problem  is  uncommon  but  serves  to  illustrate  the 
representation  of  material  property  variation  in  the  clearest  manner.  The  beam  was  analyzed  using  a  coarse 
model  of  5  elements.  Two  different  mesh  configurations  were  used  as  shown  in  Figure  4;  a  uniform  mesh 
was  adopted  to  assess  optimum  element  performance  and  a  nonuniform  mesh  was  used  to  assess  distortion 
sensitivity  to  the  simplifications  currently  incorporated  in  the  explicit  derivation.  For  simplicity,  isotropic 
material  properties  were  assumed  with  a  linear  variation  of  modulus  as  depicted  in  Figure  5. 
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Figure  4.  Cantilevered  beam  configurations:  a)  uniform  mesh,  b)  distorted  mesh. 


Table  6  depict  solutions  for  the  explicit  and  numerical  hybrid  element  formulations  together  with  results 
using  an  incompatible  displacement-based  element. 

Table  6.  Deflection  of  nonlinear  cantilevered  beam  under  end  shear  loading. 


a) Uniform  Mesh 

1  b)  Distorted  Mesh  | 

Elements 

Va 

Q8I 

Q8II 

Elements 

Va 

K88I 

PS- Explicit 

-88.7 

3587  1 

PS- Explicit 

-86.1 

KUga 

PS-Numerical 

-89.1 

PS-Numerical 

D-Based 

-89.0 

IB 

3586 

D-Based 

-89.9 

3276 

3079 

Exact 

-90.4 

4043 

3601 

Exact 

-90.4 

4043 

3601 

Comparing  the  above  results  show  an  excellent  agreement  between  the  explicit  and  numerical  hybrid  el- 
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ement  formulations  for  both  meshes;  a  small  departure  is  seen  in  the  tip  deflection  prediction  in  the  explicit 
formulation  using  the  irregular  mesh  due  to  the  incorporated  simplifications.  The  incompatible  displacement- 
based  element  demonstrates  good  results  for  both  stresses  and  displacements  using  a  uniform  mesh;  however, 
in  the  distorted  mesh,  stress  recovery  is  severely  compromised.  In  general,  the  greater  computational  cost 
of  the  displacement-based  element  makes  it  unappealing. 

REMARK  I 

An  important  aspect  of  the  above  development  involves  the  collocation  error  of  the  D  matrix  over  the 
element  domain  using  quadratic  isoparametric  shape  functions.  Because  this  matrix  is  computed  as  the 
square  root  of  the  material  stiffness  matrix,  the  individual  components  in  D  are  interpolated  to  yield  the 
approximate  relation  given  above  in  equation  (39).  A  formal  error  estimate  may  be  derived  for  arbitrary 
variation  of  material  properties  and  order  of  collocation  polynomials;  however,  a  clearer  illustration  may  be 
given  by  assuming  a  specific  variation  in  the  compliance  matrix  components  as  a  function  of  a  single  coor¬ 
dinate,  thus  simplifying  the  error  analysis  to  a  one-dimensional  demonstration.  A  linear  function  is  selected 
in  which  a  parameter,  Q,  is  used  to  set  the  magnitude  of  variation  in  the  material  component  denoted  by  Cq. 
The  resulting  square  root  distribution  is  given  by 

/(0  =  Q'^'[i  +  s(i  +  0F 

The  values  of  /(^)  at  the  evaluation  points  along  ^  6  (—1, 1)  are  depicted  in  Figure  6. 


f(5)f 


The  quadratic  isoparametric  interpolation  is  given  by 

p(0  =  +  |[(1  +  +  ^[1  -  2(1  +  0)1/2  +  (1  +  20)1/2]^2} 

An  error  measure  for  point  evaluation  may  be  given  by 

p  l/«)-i>(OI 

’ - 7W~ 

A  second  error  measure  is  associated  with  the  difference  between  the  integrated  areas  computed  by  the  exact 
solution  and  the  quadratic  collocation.  The  integral  of  the  exact  distribution  is  given  by 

m  =  j\  nm  =  - 1] 

while  the  integral  of  the  collocated  function  is  given  by 

/’(O  =  +  4(1  +  0)1/2  +  (1  +  20)1/2] 
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An  error  measure  for  the  integrcil  evaluation  may  be  given  by 


Ei  = 


m)  -  pm 
m 


Table  7  details  error  measures  at  the  points  of  maximum  expected  error  and  the  integrated  error  between 
the  two  assumed  functions. 


Table  7.  Error  measures  for  point  and  integral  evaluation  of  material  property  interpolation. 


m 

■jilJ-MW 

Ei 

■miifl 

0.0 

0.0 

0.0 

2.287E-6 

2.262E-6 

5.005E-9 

2.603E-4 

2.467E-4 

2.683E-6 

MM 

1.862E-3 

1.681E-3 

3.581E-5 

MM 

1.211E-2 

1.003E-2 

4.109E-4 

MM 

3.384E-2 

2.606E-2 

1.542E-3 

0.4 

6.739E-2 

4.871E-2 

3.710E-3 

0.5 

1.120E-1 

7.658E-2 

7.048E-3 

0.6 

1.664E-1 

1.083E-1 

1.158E-2 

0.7 

2.292E-1 

1.429E-1 

1.728E-2 

0.8 

2.993E-1 

1.794E-1 

2.406E-2 

0.9 

3.754E-1 

2.171E-1 

3.184E-2 

1.0 

4.565E-1 

2.555E-1 

4.051E-2 

The  above  table  indicates  that  the  maximum  point  collocation  error  for  a  40%  variation  in  modulus  {9  =  0.2) 
over  the  segment  is  just  over  1%.  At  ^  =  .0.5,  corresponding,  to  a  100%  variation  in  material  properties,  the 
point  error  increases  to  11.2%.  However,  even  when  the  properties  vary  by  a  factor  of  3  (^  =  1.0)  over  the 
segment,  the  integrated  error,  is  still  small,  on  the  order  of  4.05%. 

REMARK  II 

With  stresses  assumed  in  natural  coordinates,  the  procedure  described  herein  may  be  applied  directly.  For 
higher-order  elements,  the  approach  of  assuming  tensorial  stresses  with  a  contravariant  transformation  us¬ 
ing  centroidal  Jacobians  to  obtain  Cartesian  stresses  is  inaccurate  and  neccessitates  a  formulation  based 
on  stresses  assumed  a  priori  in  Cartesian  coordinates.  For  assumed  Cartesian  stress  fields,  a  fully  explicit 
derivation  may  become  overly  cumbersome  and  computationally  disadvantageous.  In  such  cases  the  basic 
methodology  is  applied  but  numerical  quadrature  of  the  sczdar  integrals  arising  in  the  weighted  inner  product 
and  computation  of  G-matrix  components  is  advocated.  In  higher-order  elements,  the  reduction  in  compu¬ 
tational  cost  afforded  by  adapting  the  present  methodology  is  expected  to  be  significantly  greater  than  that 
demonstrated  for  the  4-node  quadrilateral  element  due  to  the  larger  order  of  the  constituent  matrices.  The 
computational  savings  result  from  eliminating  the  cost  of  forming  and  inverting  the  complementary  energy 
matrix,  H,  and  by  replacing  the  numerical  quadrature  of  large-order  matrix  products  by  the  quadrature  of 
a  small  set  of  scalar  integrals. 


CONCLUSION 

The  developed  methodology  for  deriving  explicit  hybrid  element  stiffness  matrices  has  been  applied  to  the 
linear-elastic  Pian-Tong  8-node  solid  continuum  element  and  has  been  extended  to  develop  an  explicit 
formulation  for  the  nonlinear-elastic  Pian-Sumihara  quadrilateral  element.  The  3-D  element  formulation  in¬ 
corporates  a  complex  set  of  higher-order  stress  modes  yet,  through  application  of  a  simplifying  methodology, 
has  been  shown  to  permit  a  straightforward  derivation  of  explicit  algebraic  expressions  for  element  stiffness 
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coefficients.  In  comparison  with  an  underintegrated  numerical  version  of  the  element  using  a  2”^-order  Gaus¬ 
sian  quadrature  rule,  the  explicit  derivation  required  less  than  14  %  of  the  computational  cost  in  forming 
element  stiffness  matrices.  If  an  exact  integration  is  required  in  the  numerical  version,  the  explicit  derivation 
demonstrates  a  20-fold  increase  in  computational  efficiency.  The  enhancement  of  computational  efficiency  of 
hybrid  element  formulations  for  nonlinear  analysis  has  been  demonstrated  using  the  4-node  Pian-Sumihara 
quadrilatercJ  element  to  derive  explicit  formulations  acconuriodating  nonconstant  material  properties.  The 
derivation  of  an  explicit  element  stiffness  matrix  has  been  shown  to  substantially  reduce  the  computational 
cost  in  nonlinear  analysis.  The  developed  methodology  is  completely  generic  and  may  be  applied  to  any 
hybrid  element  formulation  to  reduce  the  computational  cost  in  linear  and  nonlinear  applications.  In  the 
application  to  higer-order  element  formulations,  an  assumption  of  Cartesian  stresses  lead  most  efficiently  to 
a  combination  of  numerical  evaluation  of  the  scalar  integrals  involved  in  the  orthonormalization  process  and 
in  the  integration  of  of  various  integrals  required  in  determining  components  of  the  G  matrix.  All  other  inte¬ 
grations  may  be  performed  analytically.  The  increase  in  efficiency  demonstrated  with  a  linear-order  hybrid 
element  is  expected  to  be  even  more  pronounced  in  higher-order  element  formulations  due  to  the  elimination 
of  forming  and  inverting  the  complementary  energy  matrix  and  by  replacing  the  numerical  quadrature  of 
large-order  matrix  products  with  the  analytical/numerical  integration  of  a  relatively  small  set  of  scalar  inte¬ 
grals.  The  application  of  the  above  methodology  can  be  expected  to  find  general  application  in  hybrid  and 
mixed  element  formulations  and  provide  a  significant  reduction  in  computational  cost  in  generating  element 
stiffness  matrices  for  hnear  and  nonlinear  analysis. 


APPENDIX  I 

SPECTRAL  DECOMPOSITION 


The  spectral  decomposition  of  an  orthotropic  material  compliance  matrix  is  given  by 


where  the  C  and  D  matrix  are  given  for  a  2-D  orthotropic  material  as 


’  Cii 

Cl  2 

0  ■ 

■  dn 

di2 

0 

Cl2 

C22 

0 

D  = 

di2 

d22 

0 

0 

0 

C33 

0 

0 

- 1 

The  eigenvalues  are  computed  as 

(P2  =: 

^3  = 


(~y^C22  —  2CiiC22  +  4Cp  ^  ^11  "b  ^22  +  Cii)/2 
(  >/c^2  -  2ciiC22  +  4cf2  +  cfj  +  C22  +  Cii)/2 
C33 


yielding  the  A^^^  matrix  as 


The  Q  matrix  is  defined  as 
where  the  eigenvectors  are  given  by 


_  I 

f  1 

1  .  1 

{  P2  —  C22  1 

1  1 

1  “  1 

$1  = 

1  fi-cn  . 

[  0  J 

II 

1  0  J 

II 

CO 

[cl] 

and  normalized  as 

Ni  = 
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=  Ni^i 

The  computation  of  the  3-D  distributing  matrix  using  the  symmetric  C  matrix  for  an  orthotropic  material 
is  performed  as  follows. 


’  Cll 

C12 

Cl3 

0 

0 

0 

■  du 

di2 

<^13 

0 

0 

0 

C12 

C22 

C23 

0 

0 

0 

d\2 

d22 

d23 

0 

0 

0 

Ci3 

C23 

C33 

0 

0 

0 

n  — 

di3 

<^23 

d33 

0 

0 

0 

0 

0 

0 

C44 

0 

0 

u  — 

0 

0 

0 

4/44 

0 

0 

0 

0 

0 

0 

C55 

0 

0 

0 

0 

0 

dzs 

0 

0 

0 

0 

0 

0 

C66 

0 

0 

0 

0 

0 

to 

In  the  spectral  decomposition,  eigenvalues  are  obtained  as 


<pl  =  ti  -\-t2  —  u/3 

^4 

=  C44 

(p2  =  (*-(<1  +  ^2)  —  2a/3  +  \/^(^i  —  ^2))/2 

=  C55 

<p3  —  (— (^1  +  ^2)  —  2a/3  —  V— 3(<i  —  i2))/2 

=  C66 

where 


h  = 


h 

q 

a 

b 

c 


=  (r  -  y/q^  4- 

=  ;  r  =  (a6  ^  3c)/6  -  aV27 

=  -(C33  +  C22  +  Cii) 

=  (C22  +  Cu)c33  —  C23  —  cf3  -|-  C11C22  “  ^12 

=  (C11C22  +  ^12)^33  +  C11C23  —  2C12C13C23  +  Ci3<^22 


yielding  the  matrix  as 


The  Q  matrix  is  given  by 

where  the  eigenvectors  are  given  by 


'Pa^ ^ 

Q  =  [^i|^2|^3|^4|^5|^6] 


ic22  -  Pl)iC33  -  Pl)  -  cla 

C13C23  —  C12(C33  -  <p2) 

C13C23  —  C12(C33  —  Pi) 

(Cll  -  V’2)(c33  -  P2)  -  Ci3 

C12C23  —  Ci3(C22  —  Pi) 

^  0 

^  ,  <>2  =  ^ 

C12C13  —  C23(Ci1  —  P2) 

0 

0 

0 

0  J 

0  J 

C12C23  —  Ci3(C22  — 

P3)  ' 

0  ] 

'  0  ^ 

'  0  ' 

C12C13  —  C23(C11  — 

P3) 

0 

0 

0 

(Cll  —  ¥>3)(C22  -  ^3) 
^  0 

-  />2 
—  C12 

>  ,^A=  < 

0 

C44 

^  $5  =  < 

0 

0 

>  ,  ^6  =  ^ 

0 

0 

0 

0 

C55 

0 

0 

\ 

J 

0  J 

0  J 

,  C66  , 

and  normalized  to  yield 

Ni  = 

^i  =  Nih 
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In  the  case  of  degenerate  eigenvalues,  the  associated  eigenvectors  are  discarded  and  replaced  by  vectors 
orthonormalized  to  the  independent  eigenvectors  using  the  standard  Gram-Schmidt  procedure.  In  addition, 
the  inverse  of  the  distributing  matrix  is  obtained  by  simply  replacing  the  diagonal  A  matrix  in  the  above 
with  the  following 


APPENDIX  II 

ORTHONORMAL  STRESS  MODE  COEFFICIENTS  FOR  THE  3-D 
PIAN-TONG  HEXAHEDRAL  ELEMENT 


For  3-D  coordinate  transformations,  the  determinant  of  the  Jacobian  is  given  by 


|J|  =  Si  +  S2^  +  S37/  +  S4C  +  SsT;^  +  SeC^  +  SyC??  +  +  SgT]^  +  SioC^  +  +  Si2T)^^  + 

+  sisC??^  +  +  Si7f?c^  +  +  S20»?^C^ 

where 


Si  =  <p^2 

56 

—  ¥>31  +  ¥>51  +  ¥>43 

Sll  =  2^54 

S16  =  ¥35 

S2  =  ¥>31  +  V7i2 

57 

=  ¥>23  +  9^34  +  ¥>26 

S12  =  ¥>14 

Sl7  =  ¥37 

S3  =  ¥>23  +  ¥>12 

58 

=  V>14 

Sl3  =  ¥>51 

S18  =  ¥54 

S4  =  ¥>31  "b  T23 

59 

=  ^>42 

Sl4  =  ¥>l2 

Sl9  =  ¥46 

S5  =  ¥>42  +  ¥>14  +  ¥>[2 

5l0 

=  ¥>35 

Sl5  =  ¥26 

S20  =  ¥65 

and 

p'jk  =  <ii{f>jCk  -  bkCj)  +  6,(cjajt  -  CkOj)  +  Ci{qjbk  -  0*6,).;: 

The  stress  mode  coefficients,  r,y,  are  given  by 


'■•.1  =  c,-2 

ri,2  =  —Ci2^2/^l 

n,3  =  Ci3  — 

r,_4  =  —CizXg/Xi  —  <t>\ri^2 

n,5  =  Ci5  -  (Af  r,-,3  -  (?i^r,-,i 

»'«,6  =  —4>iri,A  -  02’*i,2  —  Ci^Xg/Xx 

^1,7  =  Cji 

n-.s  =  -<^in,5  - 

n-,9  =  -  <f>2ri,4  -  <^3'*<,2  -  C,iA3/Ai 

»*»,10  =  C,-3  —  <^iri,7 

n-.ii  =  -<i>tri,s  -  <A2»**,5  -  <l>3ri,3  -  <t>\ri^i 

n.n  =  -4>\ri^g  -  (i>2rifi  -  <i)iri_4  -  4>\hi,2  -  C13A3/A1 

^*,13  —  ^16  ^2^1,7 

n,i4  =  -<?iiri.n  -  ^ir<,8  -  -  <^^r,-,3  - 

n-.is  =  -<t>\ri,\2  -  ^2''t,9  -  <t>3rifi  -  4>\ri,A  -  <i>\ri,2- 

CjsAa/Ai 

n,n  =  -<!iir<,i3  -  '^Idmo  -  «^3'’*,7 

nMS  =  -<^l»’,-,14  -  <t>2''iM  -  <t>3^ifi  -  ^4^i,5  -  <i>t'^i,3- 
n,19  =  15  -  (^2'’«,12  -  <i>3ri,9  -  <^4''t,6  “ 

<i>t^i,2  —  CuXa/Xi 


ri,28  =  Ci3 

»*i,29  =  —4>\u,2A  —  <^2’'«,20  —  <^3<^«1 

»'i,30  =  -^in,25  -  <^2'’*, 21  -  -  <A4»'.-,13- 

7.-, 31  =  -<Al7i,26  -  <^27f,22  -  <;i37.M8  “  <l>A'^i,lA- 
<^57t.ll  -  <^67i,8  -  4>7n,3  -  4>%ri,2  - 
7i,32  =  -<^l7i,27  -  <;i27t,23  “  <^37.',19  “  <^47t,15- 

12  -  <i>%n,9  -  <i>7^i,6  -  <^|7i,4  -  <^97i.2 
CjaAs/Ai 
7j,33  =  Cil 

71.34  =  — <Ai°7,_28 

71.35  =  -<fri°7,_29  -  '^2‘’7.-,24  “  <^3°7i,20  “  <i>A°Cil 

71.36  =  -'^i‘’r,-,3o  -  <A^°r,-,25  -  <;i3°r,-,2i  -  <p\°ri^i7- 

<f>5°^i,13  -  '^6°7i,10  -  'A7°7i,7 
7.-, 37  =  -^i°7,-,3i  -  4>2°ri^26  -  '^3°7,-,22  -  9i4“7,-,i8- 

<4s°7i,14  -  <^6°7,-,11  -  <4^‘’7i.8  -  «A8“7,-,5- 

7, ',38  =  -^i°7,-,32  -  <^2°7,-,27  “  4>3°ri,23  “  '?i4°7<,19- 
-  4>6°ri, 12  -  07“7i,9  -  08°7<,6- 
4>9°ri,A  -  <^1071,2  -  cnAe/Ai 
7., 39  =  C,-2 

7<,40  =  -<^P7i,33 
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1—  r,*,45  =  — 


n,27  =  “<^f^i,23 


r  I  •  r  O 

-0T»’t,23  -  <l>2n.is  -  ^in.is  -  <^4»'«M2  -  4>%ri.d- 
<i>6n,6  -  4>7‘ri,i  -  <;A|r.-,2  -  c,-4A4/Ai 

Expressions  for  the  inner  products  arising  from  orthogonalization 

6 


<^“*•.•,34  • 

<;^“»*.',36  • 

<^“^,18 
•'^i^»‘t,38  ■ 
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-  </>2^»'«,29 

-  <^“»'<,30 

-  '^2^»’»,31 

-  <^“’•<,14 

4>\ori,3  - 

-  <^2^»’i,32 

-  <^“’•<.15 
^10*’«,4  ~ 


-  <;i“n-.24  - '^“»’f,20- 


-<^“n,25 

—  <^7^ri|io  ■ 

—  ^3^»’«,26 


t  — 


^li’’i,l 


-<^4^''i,22- 


»*»,27  —  ^“»*j,23- 
^11>*«,2  —  Ci2^7/\l 


are  given  by 


where 

^!..- 
<t>h 

<f>h 

^'4  = 
<i>4 

4>\^.  - 

<i>4 

<l>\^ 

<t>4 

^kj 

4>1^a  - 

<t>\ii 

<!>{) 

<i>4 
<^1- 
<i>4 


= 

n7Ct.3(n- 

= 

n^Ci,5{ri 

= 

ngCi.siri 

= 

n9C,-,i(r,- 

= 

n^Ci^i(ri 

= 

nlci,i{ri 

= 

ngCi_3(ri 

= 

nic,-,3(r,- 

= 

nyCi^3(ri 

= 

ngCifiin 

= 

nlci^eiri 

= 

n7Cifi{ri 

n9C,M(r,- 

= 

nlci,i{ri 

= 

n^Ci^i(ri 

= 

nlci,2{ri 

= 

ngCiM^'i 

= 

■  «7C.-.2(n 

= 

«9Ct,4(r,- 

= 

n|ci.4(n 

= 

n7C.-,4(r,- 

= 

W9C.,3(»’i 

= 

nlci,3iri 

= 

n^Ci^siri 

= 

ngCi, !(»•,• 

= 

nlci,i{ri 

= 

nlci^2(ri 

=: 

n8C<.2(ri 

= 

njCi^2iri 

= 

«ioCi.3(»- 

”iic.-.6(r 

=r 

«ioC.-.6(r 

= 

«12C<,l(r 

,3-^5  +  ^»,4'^3) 

,1^5  4*  ri  2^2) 
+  n.eAs) 


,3^7  4-  n,4A4) 


,3Ai6  4-  r,-,4A7) 


c  =  E'^"..- 
»=1 


<i>4 

<i>4 

<t>4 

‘i>4 

4>4 

<i>4 

4>% 

A 

<i>4 

<i>4 

<t>4 

<i>4 


«!2<^*.2(»'i 
nfiC, 
«ioC. 

ni2C( 

2 


=  niiC, 

=  nfo< 
=  n?2C, 

=  nl^c, 

=  riloC, 
=  «12^1 
=  nfiC, 
=  nloCr 
=  n?2C, 

=  niic, 
=  «io< 
=  «13Ci 

=  n?4C, 

= 

=  Tli^Ci 

=  nl^c, 
=  njgc, 


-  «2 


nUCi 
=  nfgc, 

=  ^15^1 

=  nl^c^ 

=  TijaC, 

=  nfgC, 


».2(n 

,2(n 

*\4(n 

.4(r, 

.40 

»,3(n 

ti3(rt 

,i(n 

.lO 

.2(r, 

.2(n 

.2O 


n,7A6  4-  ri.gAr  4-  n,9A4) 

,13A6  +  ri^i4A7  4"  r,-^i5A4) 
.loAe  4"  r,',iiA7  4“  ri^i2A4) 

,7A6  +  ri^s^T  4-  n'^9A4) 
.laAe  4-  n*,HA7  4-  n^i5A4) 
,ioA6  4-  J',-,iiA7  +  r,-,i2A4) 
jAe  4-  ri,8A7  4-  rt,9A4) 


,13Ai2  +  n,i4Ai4  4-  rj-.isAs) 
,ioAi2  4*  r,-,iiAi4  4“  r,-,i2A5) 
,7Ai2  +  ri,8Ai4  4-  r,-^9A5) 
,13Ai7  4- n.nAu  4“r,*^i5A6) 
,ioAi7  4-  r, -^11  All  4-  r,’^i2A6) 
,7Ai7  4“  Ti.gAii  4*  n^gAe) 
,13Aii  4-  n^i4Ai6  4-  r,-,i5A7) 
,ioAii  4-  n^iiAie  +  ri^i2A7) 
,7Aii  4-  Ti^sAie  4-  rf^9A7) 


,2(^»,16Aio  +  n,17A6  4-  r,*,i8A7  4-  r,-,i9A4) 
,4(^*t,2oAiO  +  n^2lA6  4-  ri^22^7  4"  7*»,23A4) 
,4(n,i6Aio  4-  ri,i7A6  4-  ri^i8A7  4-  ri,i9A4) 
,3(n,24Aii  4-  ri,25Ai2  4-  ri,26Ai4  4“  rt^27A5) 
,3(^'t,2oAii  4“  r,‘,2i Ai2  +  r,*,22Ai4  4-  r,-,23A5) 
i,3(»'*,16All  4-  r,-,i7Ai2  +  r,-,i8Ai4  4-  r,-,i9A5) 
,i(rt,24Ai5  4-  r,*,25Ai7  4*  r,-,26Aii  +  ri^27^6) 
,i(»'t,2oAi5  4“  ri,2iAi7  4-  ri,22Aii  4-  ri,23A6) 
,i(^t,i6Ai5  4-  >’i,i7Ai7  4-  ^i^igAii  4-  r^^igAe) 
,2(n,24Ai3  4-  rt,25Aii  4-  ri,26Ai6  4-  r,-,27A7) 
,2(^i,2oAi3  4"  n^2lAii  4-  n,22Al6  4-  ri^23A7) 
i, 2(^*1, isAis  4-  ri^i7Aii  4-  rj-^gAis  4-  ri^i9A7) 
,i(^*,28A22  4-  ri,29Ai5  4"  r,-,3oAi7  4-  r*^3iAii4* 
n,32A6) 

’^i6^«,2(^t,28A23  4-  r,',29Ai3  4-  r,*^3oAii  4-  rj^3iAi64- 

'•*,32A7) 

^17C»,2(^i,33A21  4-  ri,34A23  4"  r,-,35Ai3  4"  r,‘,36Aii4- 
ri,37Ai6  4-  ri  3s^7) 


=  »^iiC»,i(rt,ioA6  4*  r,*,iiA7  4-  r,*,i2A4) 
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Normalizing  factors  for  the  othogonal  stress  modes  in  which  j  =  1, 2, 3  are  given  by 

=  {Ei=l(^-,3j+4A9  +  2r,-,3j4.4(r,'^3jf+5A5  -h  ri,3j+6A3)  +  ^,^,37+5^8  +  2ri,3j-{.5r,-^3y +6 ^2  + 

=  {Et=:i(^Mi+i2Aio  +  2r,‘,4y+i2(ri,4j4.i3A6  +  ri,4j-4.i4A7  +  rj,4j4.i5A4)  + 

2ri,4j  +  i3(r,',4j4.14A5  +  r,-,4j4.i5A3)  +  r?4^.^i4A8  +  2ri,4j4.i4r,'^4j4.15A2  +  ^?,4;-|.15Ai)} 

“  {Ei=l(^?, 28^18  +  2r, -^28(^,29^11  +  r,-^3oAi2  +  r,-^3iAi4  +  r,-^32A5)  +  29^10  + 

2^*1, 29(^i, 30^6  +  ^1,31^7  +  ^1,32^4)  +  ^^,30^9  +  2r,- 3o(r,',3i  As-f 
ri,32A3)  +  r?3lA8  +  2ri,3in,32A2  +  ^^32Al)}’^/^ 

—  {Ei=l(^?,33Al9  H-  2rj^33(ri^34A22  +  n‘,35Ai5  +  r,-,36Ai7  +  r,-^37Aii  +  r,-^38A6)+ 

^?,34Ai8  +  2r,'^34(ri^35Aii  +  r,',36Ai2  +  r,-,37Ai4  +  r»,38A5)  +  ^?,35Aio  + 

2r,-,35(ri^36A6  +  r,- 37A7  -f  r£^38A4)  -f  +  2 r,- ^35 (r,- 37 A 5  +  r,-^38A3)+ 

^\37A8  +  2ri,37n-,38A2  +  rf38Ai)}“^/2 

^  {E$=i(^?,39A20  +  2ri^39(r,-^4oA2l  +  r,-,4iA23  4-  r,*^42Ai3  -1-  r, -^43 All  +  J^t,44Ai6  +  ?’i,45A7)+ 

^i\4oAi9  +  2r,-^4o(r,-^4iA22  +  ri^42Ai5  +  7*»,43Ai7  +  ri,44Aii  4-  r,-^45A6)  4-  r?4iAi8H- 

2r,-,4i(r,-^42Aii  4-  r't,43Ai2  4*  ^t,44Ai4  +  ^i,45A5)  4"  ^,^,42Aio  4-  2r,-^42(r,-^43A64- 
n,44A7  4-  ri,45A4)  4-  ^1^,43 Ag  +  2rj^43(r,-^44 A5  4-  r, •, 45^34“ 

^»\44A8  4-  2r,-^44r,-,45A2  4- 

The  integrals  arising  in  the  weighted  inner  product  evaluate  to 


Ai 

=: 

8(3si  4"  58  4-  59  4-  5io)/3 

Ai3 

= 

8(15s2  4"  5si4  4*  9si8)/135 

A2 

— 

8(3s2  4“  5i4  +  5i6)/9 

Ai4 

= 

8(15s3  4*  9si2  4-  5si7)/135 

A3 

= 

8(3s3  4- 5i2  4- 5i7)/9 

Ai5 

= 

8(15s3  4“  5^12  4"  9si7)/135 

A4 

= 

8(3s4  4"  5i3  4-  Si5)/9 

Ai6 

= 

8(15s4  4"  9si3  4"  5si5)/135 

A5 

= 

8(3s5  4"  52o)/27 

Ai7 

8(15s4  -b  5si3  4-  9si5)/135 

Ag 

8(3s7  4-  5i8)/27 

Ai8 

= 

8(15si  4"  9s8  4“  9s9  4-  5sio)/135 

A7 

=: 

8(3s6  4"  5i9)/27 

Ai9 

= 

8(15si  4“  5^8  4*  9s9  4"  9sio)/135 

As 

= 

8(15si  4"  9^8  4*  5s9  4-  5sio)/45 

A20 

= 

8(15si  4-  9^8  4*  Ssg  4-  9sio)/135 

A9 

8(15si  4"  Sss  4-  9s9  4"  5sio)/45 

A21 

= 

8(5s5  +  3s2o)/135 

Aio 

= 

8(15si  4“  5s8  4-  5s9  4-  9sio)/45 

A22 

= 

8(5s6  +  3si9)/135 

All 

Ai2 

= 

8sii/27 

8(15s2  4“  95i4  4*  5si6)/135 

A23 

— 

8(5s7  4-3si8)/135 

APPENDIX  III 

ORTHONORMAL  STRESS  MODE  COEFICIENTS  FOR  THE 
NONLINEAR  PIAN-SUMIHARA  ELEMENT 


ni 

«i+6 

«i+9 

’^i+12 

ni6 

ni7 

^18 


The  stress  mode  coefficients,  pj'y,  are  given  by 


in  which 


Pl3 

= 

ni 

P42  = 

n4Cii 

P43  = 

— n4CiiA2/Ai 

P26 

= 

ni 

P45  = 

n4C21 

P46  = 

— n4C2iA2/Ai 

P39 

=z 

ni 

P48  = 

”4^31 

P49  = 

— n4C3iA2/Ai 

Psi 

=r 

n5Ci2 

P52  = 

-n5Cii<^ 

P53  = 

— ris^i 

P54 

”5^22 

P55  = 

-n5C2l4> 

P56  = 

-ns  ^2 

P57 

— 

”5^32 

Pss  = 

—  «5C31<A 

P59  = 

— ns^s 

<j)  ==  — n4(ci2Cii  +  C22C21  +  C32C3i)A2A3/Ai 

9i  =  (C12A3  H- (^ciiA2)/Ai 

$2  =r  (C22A3  4- <^C2iA2)/Ai 
^3  =  (^^32A3  4- <^C3iA2)/Ai 


1/2 


20 


p 


and  normalizing  factors  are  given  by 

«.  =  A.-*'’ 


7X4  —  [(<^11  +  C21  -h  C3i)(A5  —  A2^/Ai)] 

3 

7x5  =  [y^(c^^2  2c,*20iA3  +  +  2(l>Cii0i\2  + 


izzl 

The  determinant  of  the  Jacobian  for  2-D  transformations  is  given  by 

1J|  =  </o  +  Ji^  + 

where 

Jo  =  01^2  “  \  Jl  ~  ^  <13^1  t  J2  —  ^3^2  “  0>2^z 

The  scalar  integrals  arising  in  the  inner  product  to  evaluate  to 


Ai  =  J  J  ^ \J\d^dT,  =  4Jo  =  y  ^  j  [\3\^rfd^dT)  =  0 
A2  =  j'j'^\Jmdij  =  4Ji/3  A5  =  j'j'[\3\ed^dT,  =  Uo/2 
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